Técnicas en aprendizaje estadístico - Introduction to Statistical Learning - ISLR
#Librerías
library(tidyverse) # Manipulación, limpieza y gráficos
library(ISLR) # Bases de datos
library(plotly) # Interactividad gráfica
library(ggthemes) # Presentación
library(GGally) # Correlación
library(yardstick) # Métricas
library(MASS) # Clasificadores
library(class) # Clasificadores
library(tree) #CART - ISLR
library(broom) # Resultados tidy de modelos
library(latex2exp) # Expresiones LaTex en las gráficasSección 4.7 - Ejercicio 10
a) Análisis descriptivo
El conjunto de datos Weekly contiene porcentajes de retorno semanales del índice S&P entre los años 1990 y 2010 con 1089 observaciones.
Year: es el año en el que fue obtenida la observación.Lag 1-5: Retornos porcentuales de las 5 semanas anteriores, respectivamente.Volume: Volumen de acciones intercambiadas (Número promedio de acciones diarias intercambiadas en billones).Today: Retorno porcentual de la semana.Direction: Indica si el mercado tuvo un retorno positivo o negativo en esta semana.
Encabezados y primeras 6 observaciones:
## # A tibble: 6 x 9
## Year Lag1 Lag2 Lag3 Lag4 Lag5 Volume Today Direction
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 1990 0.816 1.57 -3.94 -0.229 -3.48 0.155 -0.27 Down
## 2 1990 -0.27 0.816 1.57 -3.94 -0.229 0.149 -2.58 Down
## 3 1990 -2.58 -0.27 0.816 1.57 -3.94 0.160 3.51 Up
## 4 1990 3.51 -2.58 -0.27 0.816 1.57 0.162 0.712 Up
## 5 1990 0.712 3.51 -2.58 -0.27 0.816 0.154 1.18 Up
## 6 1990 1.18 0.712 3.51 -2.58 -0.27 0.154 -1.37 Down
Volumen promedio de acciones diarias intercambiadas en el tiempo
Se evidencia un crecimiento constante con un alza grande en el año 2005 y una estabilización de la trayectoria alrededor del año 2007. Se añadió un ruido a la posición de cada punto para poder evidenciar la densidad, además de transparencia. Se utilizó el método mgcv::gam para ajustar la tendencia.
p <- weekly %>%
ggplot(aes(x = Year, y = Volume)) +
geom_jitter(alpha = 0.6) +
geom_smooth() +
labs(x = NULL,
y = "Volumen (Billones)") +
theme_economist()
p_box <- weekly %>%
ggplot(aes(x = as_factor(Year), y = Volume, color = Today)) +
geom_boxplot() +
theme_economist() +
theme(axis.text.x = element_text(angle = 45)) +
labs(x = NULL,
y = "Volumen (Billones)")
ggplotly(p, width = 800, height = 500)El incremento en la volatilidad del índice acompañó su crecimiento a partir del 2007.
Correlación entre las variables
No se encontraron relaciones lineales claras entre las variables. Se presentan comportamientos esperados de la variable Today con la variable categórica Direction representada en el color, donde los valores negativos están asociados a un decrecimiento.
b) ¿Podemos predecir si el índice subirá o no mediante regresión logística?
Utilicemos regresión logística para predecir si habrá un retorno positivo del mercado en esta semana. Para ello, se utilizará el retorno porcentual de las 5 últimas semanas y el volumen sin separar un conjunto de entrenamiento y de prueba (Modelo sobreajustado).
market_logistic <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = weekly, family = binomial)
summary(market_logistic)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = binomial, data = weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.26686 0.08593 3.106 0.0019 **
## Lag1 -0.04127 0.02641 -1.563 0.1181
## Lag2 0.05844 0.02686 2.175 0.0296 *
## Lag3 -0.01606 0.02666 -0.602 0.5469
## Lag4 -0.02779 0.02646 -1.050 0.2937
## Lag5 -0.01447 0.02638 -0.549 0.5833
## Volume -0.02274 0.03690 -0.616 0.5377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1486.4 on 1082 degrees of freedom
## AIC: 1500.4
##
## Number of Fisher Scoring iterations: 4
Si realizamos una prueba de hipótesis formal para la significancia de los parámetros, encontraremos que la probabilidad de rechazar la hipótesis nula es muy grande para todas las variables predictoras excepto para el Lag2 donde al parecer hay significancia. El valor en el tiempo \(t-2\) parece tener una relación con el tiempo \(t\).
c) Matriz de confusión de la regresión logística y calidad de la predicción del modelo sobreajustado
## Up
## Down 0
## Up 1
market_logistic_class <- tibble(value = predict(market_logistic, type = "response"),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly$Direction)
market_logistic_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.561
## 2 kap binary 0.0350
Se obtiene una tasa de clasificación correcta accuracy del 56.1% utilizando regresión logística. Es importante interpretar este resultado correctamente, ya que la predicción obtenida indica que el mercado subió 557 días y bajó 54 días, y que no se está prediciendo el valor en el tiempo \(t+1\).
Cada predicción individual se debe interpretar de la siguiente manera: “Dadas las 5 semanas anteriores y el volumen de ésta, el mercado tendrá un retorno positivo”. La mayor cantidad de errores se cometió prediciendo que el mercado iba a subir y en realidad bajó, esta situación ocurrió 430 veces.
d) Modelo de regresión logística sólo con Lag2 y conjunto de prueba
weekly_train <- weekly %>%
filter(Year <= 2008)
weekly_test <- weekly %>%
filter(Year > 2008)
lag2_logistic <- glm(Direction ~ Lag2, data = weekly_train, family = "binomial")
summary(lag2_logistic)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = weekly_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4
lag2_class <- tibble(value = predict(lag2_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
lag2_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.625
## 2 kap binary 0.141
Al igual que el modelo anterior, se tienen fallas al intentar predecir las bajas del mercado, sin embargo, este se validó en un conjunto de prueba de 104 observaciones desconocidas para el modelo. El accuracy en este caso es del 62.5%. Los resultados son sorpresivamente mejores que en el modelo sobreajustado del literal b).
d) Modelo LDA (Linear Discriminant Analysis) para predecir la tendencia del mercado
## Call:
## lda(Direction ~ Lag2, data = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162
lag2_lda_class <- tibble(pred = predict(lag2_lda, newdata = weekly_test)$class,
real = weekly_test$Direction)
lag2_lda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.625
## 2 kap binary 0.141
Se obtienen las mismas clasificaciones para el modelo LDA y el modelo de regresión logística. Como el coeficiente discriminante lineal es Lag2 = 0.4414, es de esperar que en ambas categorías la gráfica de éstos sea similar. La matriz de confusión es igual a la del modelo logístico.
f) Modelo QDA (Quadratic Discriminant Analysis) para predecir la tendencia del mercado
## Call:
## qda(Direction ~ Lag2, data = weekly_train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
lag2_qda_class <- tibble(pred = predict(lag2_qda, newdata = weekly_test)$class,
real = weekly_test$Direction)
lag2_qda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo QDA con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.587
## 2 kap binary 0
El clasificador QDA clasificó todos los valores como “Up”, si bien es una estrategia que otorga resultados buenos en términos porcentuales con respecto al accuracy, se espera que el modelo proponga en algunos casos predicciones de baja del mercado.
g) Modelo KNN (K-Nearest Neighbors) para predecir la tendencia del mercado
set.seed(1)
lag2_knn <- knn(train = as.matrix(weekly_train$Lag2), test = as.matrix(weekly_test$Lag2), cl = weekly_train$Direction, k = 3)
lag2_knn_class <- tibble(pred = lag2_knn,
real = weekly_test$Direction)
lag2_knn_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN con Lag2 como único predictor")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.548
## 2 kap binary 0.0453
El clasificador \(KNN\) tiene una tasa de clasificación correcta del 50% utilizando \(k = 1\), sin embargo, al utilizar \(k = 3\) se obtiene un 54.8%. Este clasificador no tuvo inconvenientes en proponer predicciones variadas según el valor anterior del mercado.
h) ¿Qué métodos producen los mejores resultados para predecir la tendencia del mercado utilizando Lag2?
Los métodos de regresión logística y análisis de discriminante lineal obtuvieron los mejores resultados, con un accuracy del 62.5%. Se decide favorecer al modelo de regresión logística dada la facilidad de interpretación de los coeficientes y las facilidades gráficas de comunicación que permite la curva logística.
g <- augment_columns(lag2_logistic, type.predict = "response", newdata = weekly_test) %>%
arrange(.fitted) %>%
rownames_to_column() %>%
mutate(rowname = as.integer(rowname)) %>%
inner_join(lag2_class, by = c(.fitted = "value")) %>%
ggplot(aes(x = rowname, y= .fitted, color = Direction, label = Year, label1 = pred)) +
geom_point(alpha = 0.5, shape = 1, stroke = 2) +
geom_hline(yintercept = 0.5) +
theme_light() +
xlab("Índice") +
ylab("Probabilidad de que el mercado suba")
ggplotly(g, width = 800, height = 500)i) Mejora del mejor modelo y variables adicionales de interés
Selección Step-wise con la función stepAIC sin interacciones
full_logistic <- glm(Direction ~ .-Year-Today, data = weekly_train, family = "binomial")
stepAIC(full_logistic, direction = "both", trace = FALSE)##
## Call: glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = weekly_train)
##
## Coefficients:
## (Intercept) Lag1 Lag2
## 0.21109 -0.05421 0.05384
##
## Degrees of Freedom: 984 Total (i.e. Null); 982 Residual
## Null Deviance: 1355
## Residual Deviance: 1347 AIC: 1353
step_logistic <- glm(Direction ~ Lag1 + Lag2, data = weekly_train, family = "binomial")
step_class <- tibble(value = predict(step_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
step_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.577
## 2 kap binary 0.0350
El modelo seleccionado por Stepwise selection tiene un accuracy menor que el de un solo predictor, sin embargo, su AIC es el mismo. El AIC evalúa únicamente el ajuste, por lo que se recomendaría el modelo de una sola variable predictora. Es importante resaltar que el modelo de menor AIC podría ser mejor en un conjunto de prueba diferente.
Mejor modelo por tanteo de interacciones
Se proponen variables con interacciones y se selecciona el modelo con mayor accuracy en el conjunto de prueba, en este caso utilizando la variable Lag2 y la interacción entre Lag1, Lag 4 y Volume
final_logistic <- glm(Direction ~ Lag2 + Lag1:Lag4:Volume, data = weekly_train, family = "binomial")
summary(final_logistic)##
## Call:
## glm(formula = Direction ~ Lag2 + Lag1:Lag4:Volume, family = "binomial",
## data = weekly_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.493 -1.265 1.023 1.089 1.446
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2047966 0.0643540 3.182 0.00146 **
## Lag2 0.0553928 0.0291831 1.898 0.05768 .
## Lag1:Lag4:Volume 0.0007316 0.0014747 0.496 0.61982
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.3 on 982 degrees of freedom
## AIC: 1356.3
##
## Number of Fisher Scoring iterations: 4
final_class <- tibble(value = predict(final_logistic, type = "response", newdata = weekly_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "Up",
TRUE ~ "Down")), "Down", "Up"),
real = weekly_test$Direction)
final_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.644
## 2 kap binary 0.179
Con base en los resultados anteriores se selecciona el modelo por tanteo por su accuracy superior, sin embargo, esto puede ser un sobreajuste del conjunto de prueba; por ello, si este modelo fuera a producción se recomienda realizar validación cruzada para identificar si la interacción de variables añadida es significativa. Si no es lo es, el modelo seleccionado por Stepwise selection se considera apropiado.
Sección 4.7 - Ejercicio 11
En este apartado, desarrollaremos un modelo para predecir cuando un carro tiene alto o bajo rendimiento de combustible respecto al kilometraje.
Inicialmente observemos las características de este conjunto de datos. El conjunto de datos Auto contiene información sobre 392 vehículos cuyas variables son:
mpg: millas por galón.cylinders: numero de cilindros (entre 4 y 8)displacement: desplazamiento del motor en pies.horsepower: caballos de fuerza del motor.weight: peso del vehículo en libras.acceleration: tiempo que tarda en acelerar de 0 a 60 medido en segundos.year: modelo del auto (módulo 100)origin: origen del vehiculo.name: nombre del vehículo.
Se mostrarán los primeros 6 datos junto con su encabezado.
## # A tibble: 6 x 9
## mpg cylinders displacement horsepower weight acceleration year origin
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 18 8 307 130 3504 12 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11 70 1
## 4 16 8 304 150 3433 12 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10 70 1
## # ... with 1 more variable: name <fct>
a) Creación de una variable respuesta binaria
A continuación crearemos la variabel mpg01 que tomará el valor de 1 si la variable mpg esta por encima de la mediana, y 0 en caso contrario.
median_mpg <- median(auto$mpg)
auto <- auto %>% mutate(
mpg01 = case_when(
mpg >= median_mpg ~ 1,
mpg < median_mpg ~ 0
),
mpg01 = as.factor(mpg01),
origin = as.factor(origin)
)
# Retiramos la variable mpg
auto <- auto %>% dplyr::select(-mpg)
head(auto)## # A tibble: 6 x 9
## cylinders displacement horsepower weight acceleration year origin name
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <fct>
## 1 8 307 130 3504 12 70 1 chev~
## 2 8 350 165 3693 11.5 70 1 buic~
## 3 8 318 150 3436 11 70 1 plym~
## 4 8 304 150 3433 12 70 1 amc ~
## 5 8 302 140 3449 10.5 70 1 ford~
## 6 8 429 198 4341 10 70 1 ford~
## # ... with 1 more variable: mpg01 <fct>
b) Relaciones entre las variables
A continuación realicemos un análisis descriptivo para observar el comportamiento de las variables regresoras respecto a la variable mpg01
ggpairs(data = auto, columns = 1:6, mapping = aes(color = mpg01), legend = 1, lower = list(continuous = wrap("points", alpha = 0.6))) +
theme_light() +
theme(legend.position = "bottom")Por otro lado, demos una visualización para las variable origin, la cual es categórica.
ggpairs(data = auto, columns = c(7, 9), mapping = aes(color = mpg01), legend = 1, lower = list(continuous = wrap("points", alpha = 0.6))) +
theme_light() +
theme(legend.position = "bottom")Tomando en cuenta las gráficas, podemos afirmar que las variables cylinders, displacement, horsepower y weight son buenas variables cuantitativas para realizar una clasificación de la variable mpg01, pues al graficar discriminando por clases se puede observar que hay una particion notoria entre los carros de alto consumo y bajo consumo respecto a estas 4 variables. Note también que la variabele origen puede ser una buena variable predictora para mpg01 cuando el tipo de carro es americano.
c) División de conjuntos de entrenamiento y validación
A continuación dividiremos el conjunto Auto en dos conjuntos, uno de validación y otro de entrenamiento dejando el 75% de los datos para entrenar los modelos.
d) Modelo LDA
A continuación realizaremos un modelo LDA considerando las varables que fueron extraidas en el análisis descriptivo y que se consideraron importantes a la hora de realizar una clasificación para la variable respuesta mpg01.
lda_model = lda(mpg01 ~ cylinders + displacement + horsepower + weight + origin, data = auto_train)
lda_model## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight +
## origin, data = auto_train)
##
## Prior probabilities of groups:
## 0 1
## 0.5204082 0.4795918
##
## Group means:
## cylinders displacement horsepower weight origin2 origin3
## 0 6.745098 270.0784 128.03922 3595.412 0.05882353 0.04575163
## 1 4.163121 115.5887 79.26241 2339.532 0.26241135 0.38297872
##
## Coefficients of linear discriminants:
## LD1
## cylinders -0.4317395853
## displacement -0.0018236085
## horsepower 0.0018307395
## weight -0.0007360167
## origin2 0.3722004164
## origin3 0.5312690065
Grafiquemos a continuación el resultado del modelo.
Ahora construyamos la matriz de confusión para el modelo con el conjunto de datos de validación para observar el desempeño.
mpg_lda_class <- tibble(pred = predict(lda_model, auto_test %>% dplyr::select(-mpg01))$class ,
real = auto_test$mpg01)
mpg_lda_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.867
## 2 kap binary 0.729
Note entonces que la tasa de clasificación correcta es de un 86.73%, lo que muestra que el modelo tiene una alta calidad.
e) Modelo QDA
Con el fin de realizar una comparativa con el modelo anterior, realicemos un modelo QDA para predecir la varibale mpg01 usando las variables más significativas encontradas en el análisis descriptivo.
qda_model = qda(mpg01 ~ cylinders + displacement + horsepower + weight + origin, data = auto_train)
qda_model## Call:
## qda(mpg01 ~ cylinders + displacement + horsepower + weight +
## origin, data = auto_train)
##
## Prior probabilities of groups:
## 0 1
## 0.5204082 0.4795918
##
## Group means:
## cylinders displacement horsepower weight origin2 origin3
## 0 6.745098 270.0784 128.03922 3595.412 0.05882353 0.04575163
## 1 4.163121 115.5887 79.26241 2339.532 0.26241135 0.38297872
Ahora construyamos la matriz de confusión para el modelo QDA con el conjunto de datos de validación.
mpg_qda_class <- tibble(pred = predict(qda_model, auto_test %>% dplyr::select(-mpg01))$class ,
real = auto_test$mpg01)
mpg_qda_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo QDA")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.867
## 2 kap binary 0.729
Note entonces que la tasa de clasificación correcta es de un 86.73%, y en comparación con el modelo LDA, se puede decir que tienen un desempeño muy parecido de acuerdo a la metrica de la tasa de clasificación correcta.
f) Modelo de regresión logística
Por otra parte, construiremos un modelo de regresión logística para realizar una predicción de la variable mpg01.
logistic_model <- glm(mpg01 ~ cylinders + displacement + horsepower + weight + origin, data = auto_train, family = binomial)
summary(logistic_model)##
## Call:
## glm(formula = mpg01 ~ cylinders + displacement + horsepower +
## weight + origin, family = binomial, data = auto_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.53196 -0.25902 -0.00542 0.33608 3.03362
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.9583896 2.1576414 5.542 2.98e-08 ***
## cylinders -0.1003094 0.4308485 -0.233 0.81590
## displacement -0.0117784 0.0122590 -0.961 0.33665
## horsepower -0.0546728 0.0176520 -3.097 0.00195 **
## weight -0.0015796 0.0009063 -1.743 0.08135 .
## origin2 0.1646890 0.6822990 0.241 0.80927
## origin3 0.6067756 0.6915181 0.877 0.38024
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 407.08 on 293 degrees of freedom
## Residual deviance: 152.08 on 287 degrees of freedom
## AIC: 166.08
##
## Number of Fisher Scoring iterations: 7
Calculemos de igual manera la matriz de confusión para realizar una comparación con los métodos previos usando los datos de entrenamiento.
predictions_prob_glm <- predict(logistic_model, auto_test %>% dplyr::select(-mpg01), type = "response")
predictions_glm <- rep(0, nrow(auto_test))
predictions_glm[predictions_prob_glm > 0.5] <- 1
mpg_logistic_class <- tibble(pred = as.factor(predictions_glm),
real = auto_test$mpg01)
mpg_logistic_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo de regresión logística")test_results_logistic <- mpg_logistic_class %>%
metrics(truth = real, estimate = pred)
test_results_logistic## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.888
## 2 kap binary 0.772
De acuerdo a este resultado, se puede ver que todos los métodos tienen un puntaje muy similar entre ellos pues en este caso, se puede observar una tasa de clasificación correcta del 88.78%.
g) Modelo \(KNN\)
Realicemos a continuación un modelo de \(KNN\) para predecir la variable mpg, en este caso, debemos probar con varios casos de \(K\) hallar un óptimo. Para esto, realizaremos una función que reciba como parámetro un numero \(k\) y esta función retornará la tasa de clasificación correcta para el \(k\) dado, realizando la clasificación sobre el conjunto de prueba.
get_accuracy_knn <- function(k) {
set.seed(10)
knn_model <- knn(auto_train %>% dplyr::select(-mpg01, -name),
auto_test %>% dplyr::select(-mpg01, -name),
auto_train$mpg01,
k)
knn_class <- tibble(pred = as.factor(knn_model), real = auto_test$mpg01)
results_knn <- knn_class %>%
metrics(truth = real, estimate = pred)
return(results_knn$.estimate[1])
}Ahora tomaremos varios valores de \(k\) en un rango entre 1 y 60, y para estos valores, computaremos su tasa de clasificación correcta y los graficaremos para observar cual es el mejor valor de \(k\) en este rango.
k_values = 1:60
accuracy_results <- tibble(k = k_values, accuracy = sapply(k_values, get_accuracy_knn))
g <- ggplot(data = accuracy_results) +
geom_line(mapping = aes(x = k, y = accuracy)) +
xlab("K") +
ylab("Tasa de clasificación correcta") +
theme_light()
ggplotly(g, width = 800, height = 500)Una vez observada la gráfica, encontraremos el mejor valor de \(k\) que mejora la tasa de clasificación correcta, y además para este valor de \(k\) realizaremos la matriz de confusión para ver el desempeño del modelo de manera global.
set.seed(10)
best_k <- min(
accuracy_results[which(accuracy_results$accuracy == max(accuracy_results$accuracy)), "k"]
)
knn_best_model <- knn(auto_train %>% dplyr::select(-mpg01, -name),
auto_test %>% dplyr::select(-mpg01, -name),
auto_train$mpg01,
best_k)
mpg_knn_class <- tibble(pred = as.factor(knn_best_model), real = auto_test$mpg01)
test_results_knn <- mpg_knn_class %>%
metrics(truth = real, estimate = pred)
mpg_knn_class %>% conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN")Concluimos entonces que para los valores de \(k\) probados, el mejor accuracy obtenido fue de 89.8%, el cual sigue siendo muy similar a los puntajes anteriores y desde un punto de vista global, estos puntajes denotan que los modelos tienen una muy buena calidad.
Sección 4.7 - Ejercicio 12
a) Implementación de la función Power()
A continuación implementaremos una función Power() que calculará la potencia \(2^3\)
b) Implementación de la función Power2() generalizada
Ahora, implementaremos la función Power2() que recibe dos parámetros x y a y computará \(x^a\)
c) Usos de la funcion Power2()
En el siguiente bloque de código calcularemos \(10^3\).
## [1] 1000
Ahora calcularemos el valor de \(8^{17}\).
## [1] 2.2518e+15
Y finalmente calcularemos el valor de \(131^3\).
## [1] 2248091
d) Implementación de la función Power3()
Ahora implementaremos una versión mucho más general para realizar las potencias que esta basada en la función Power2() pero en este caso, la función no imprimirá el resultado sino que retornara dicho valor como un objeto de R
e) Usando Plot3() para graficar la función \(f(x) = x ^ 2\)
En este bloque de código, graficaremos \(f(x) = x ^ 2\) usando una escala logaritmica para el eje \(y\)
ggplot(data = tibble(x = 1:10, y = sapply(1:10, Power3, 2))) +
geom_line(mapping = aes(x = x, y = y)) +
xlab("x") +
ylab(TeX("x^2")) +
scale_y_log10()f) Implementación de una forma generalizada para graficación
A continuación implementaremos un código para graficar la función \(f(x) = x ^ a\) tal que \(x \in I\) para un intervalo \(I \subseteq \mathbb{R}\).
PlotPower <- function(interval, a) {
ggplot(data = tibble(x = interval, y = sapply(interval, Power3, a))) +
geom_line(mapping = aes(x = x, y = y)) +
xlab("x") +
ylab(TeX(paste0("x^", toString(a)))) +
labs(title = "Gráfica de la función y = f(c)")
}
PlotPower(-10:10, 7)Sección 4.7.1 - Ejercicio 13
Utilizaremos el conjunto de datos Boston para predecir si un suburbio tiene un porcentaje de crimen mayor o menor que la media. El conjunto de datos tiene las siguientes variables:
crim: Ratio de crimen per cápita por pueblo.zn: Proporción de zonas residenciales asignadas para lotes mayores a 25.000 sq.ft.indus: Proporción de negocios no comerciales en acres por pueblo.chas: 1 si está en trayecto del río Charles y 0 en otro caso.nox: Concentración de óxidos de nitrógeno (partes por millón).rm: Número promedio de cuartos por hogar.age: Proporción de unidades ocupadas construidas antes de 1940.dis: Media ponderada de las distancias a 5 de los centros de empleo de Boston.rad: Índice de accesibilidad a autopistas periféricas.tax: Valor total de los impuestos por cada $10.000ptratio: Ratio de pupilos-profesores por pueblo.black: \(1000(Bk - 0.63)^2\) dónde \(Bk\) es la proporción de negros por pueblo.lstat: Menor estatus de la población (porcentaje).medv: Valor medio de una propiedad ocupada en $1000s.
Creación de la variable de interés
boston <- as_tibble(Boston) %>%
mutate(crimen = factor(case_when(crim > median(crim) ~ "alto",
TRUE ~ "bajo")))La mediana de la variable crim es 0.25651, se creará la variable crimen con niveles “alto” si es mayor a la mediana y “bajo” en otro caso.
Análisis descriptivo
Las variables de edad e indus presentan comportamientos discriminantes de interés, en la variable edad, los valores mayores a 70 tienen una gran prevalencia en la categoría alto y menores a este valor en la categoría bajo, por lo que al incluirla en uno de los modelos de clasificación se espera que sea significativa.
Este mismo comportamiento se identifica en la variable indus, donde los valores menores a 15 parecen tener una gran prevalencia en la categoría bajo, y entre 18 y 22 se encuentra la mayor cantidad de alto.
crim_hist <- boston %>%
ggplot(aes(x = age, fill = crimen)) +
geom_histogram(color = "black") +
labs(title = "Proporción de unidades ocupadas construidas antes de 1940.")
ggplotly(crim_hist)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
indus_hist <- boston %>%
ggplot(aes(x = indus, fill = crimen)) +
geom_histogram(color = "black") +
labs(title = "Proporción de negocios no comerciales en acres por pueblo.")
ggplotly(indus_hist)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Regresión logística
Construyamos inicialmente un clasificador únicamente con las dos características identificadas en el análisis descriptivo para evaluar las suposiciones. Realicemos también la separación en los conjutos de prueba y de entrenamiento.
set.seed(1995)
boston_train <- boston %>% sample_frac(size = 0.8)
boston_test <- boston %>% anti_join(boston_train)## Joining, by = c("crim", "zn", "indus", "chas", "nox", "rm", "age", "dis", "rad", "tax", "ptratio", "black", "lstat", "medv", "crimen")
logistic_boston_manual <- glm(crimen ~ age + indus, data = boston_train, family = "binomial")
summary(logistic_boston_manual)##
## Call:
## glm(formula = crimen ~ age + indus, family = "binomial", data = boston_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6774 -0.5407 -0.3414 0.5724 2.6981
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.727028 0.491097 9.625 < 2e-16 ***
## age -0.047380 0.007062 -6.709 1.96e-11 ***
## indus -0.131905 0.025355 -5.202 1.97e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 561.25 on 404 degrees of freedom
## Residual deviance: 337.58 on 402 degrees of freedom
## AIC: 343.58
##
## Number of Fisher Scoring iterations: 5
Efectivamente ambos parámetros son altamente significativos. evaluemos la calidad de las predicciones.
boston_manual_class <- tibble(value = predict(logistic_boston_manual, type = "response", newdata = boston_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "bajo",
TRUE ~ "alto")), "alto", "bajo"),
real = boston_test$crimen)
boston_manual_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con age y indus")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.772
## 2 kap binary 0.547
Los resultados son aceptables, sin embargo, evaluemos el modelo step-wise propuesto para identificar otras variables de interés.
boston_full_model <- glm(crimen ~ . -crimen - crim, data = boston_train, family = "binomial")
summary(boston_full_model)##
## Call:
## glm(formula = crimen ~ . - crimen - crim, family = "binomial",
## data = boston_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.4328 -0.0032 0.0000 0.1377 1.6830
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 27.787319 8.586662 3.236 0.00121 **
## zn 0.062198 0.033850 1.837 0.06614 .
## indus 0.068162 0.055431 1.230 0.21882
## chas -1.101626 0.829059 -1.329 0.18393
## nox -46.572894 8.258444 -5.639 1.71e-08 ***
## rm 0.079837 0.774442 0.103 0.91789
## age -0.033585 0.014415 -2.330 0.01981 *
## dis -0.659126 0.244784 -2.693 0.00709 **
## rad -0.663217 0.167363 -3.963 7.41e-05 ***
## tax 0.006860 0.003033 2.262 0.02372 *
## ptratio -0.458391 0.151742 -3.021 0.00252 **
## black 0.035793 0.014010 2.555 0.01062 *
## lstat -0.006730 0.058173 -0.116 0.90790
## medv -0.141170 0.073683 -1.916 0.05538 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 561.25 on 404 degrees of freedom
## Residual deviance: 160.46 on 391 degrees of freedom
## AIC: 188.46
##
## Number of Fisher Scoring iterations: 9
##
## Call: glm(formula = crimen ~ zn + nox + age + dis + rad + tax + ptratio +
## black + medv, family = "binomial", data = boston_train)
##
## Coefficients:
## (Intercept) zn nox age dis
## 26.279813 0.070380 -41.630025 -0.033514 -0.621363
## rad tax ptratio black medv
## -0.762505 0.008532 -0.402343 0.031772 -0.134618
##
## Degrees of Freedom: 404 Total (i.e. Null); 395 Residual
## Null Deviance: 561.2
## Residual Deviance: 163.3 AIC: 183.3
boston_stepwise <- glm(formula = crimen ~ zn + nox + age + dis + rad + tax + ptratio + black + medv,
family = "binomial",
data = boston_train)
boston_stepwise_class <- tibble(value = predict(boston_stepwise, type = "response", newdata = boston_test),
pred = fct_relevel(as_factor(case_when(value > 0.5 ~ "bajo",
TRUE ~ "alto")), "alto", "bajo"),
real = boston_test$crimen)
boston_stepwise_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para la regresión logística con el modelo stepwise")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.901
## 2 kap binary 0.800
boston_plot <- augment_columns(boston_stepwise, type.predict = "response", newdata = boston_test) %>%
arrange(.fitted) %>%
rownames_to_column() %>%
mutate(rowname = as.integer(rowname)) %>%
inner_join(boston_stepwise_class, by = c(.fitted = "value")) %>%
ggplot(aes(x = rowname, y= .fitted, color = crimen, label1 = pred)) +
geom_point(alpha = 0.5, shape = 1, stroke = 2) +
geom_hline(yintercept = 0.5) +
theme_light() +
xlab("Índice") +
ylab("")
ggplotly(boston_plot, width = 800, height = 500)El modelo stepwise tiene una mejora considerable en las métricas al incluir otras variables adicionales, sin embargo, se destaca que el poder discrimatorio de las variables adicionales con respecto al modelo inicial propuesto sólo mejoró la predicción en un 12.9%. En un ambiente de producción se recomendaría evaluar el costo de utilizar estas variables en el modelo versus el beneficio de predecir correctamente.
Modelo LDA
boston_lda <- lda(crimen ~ zn + nox + age + dis + rad + tax + ptratio + black + medv, data = boston_train)
boston_lda## Call:
## lda(crimen ~ zn + nox + age + dis + rad + tax + ptratio + black +
## medv, data = boston_train)
##
## Prior probabilities of groups:
## alto bajo
## 0.5111111 0.4888889
##
## Group means:
## zn nox age dis rad tax ptratio
## alto 1.178744 0.6392367 86.25942 2.528803 14.787440 507.3575 18.97874
## bajo 22.209596 0.4703273 51.19040 5.149906 4.136364 306.2677 17.84192
## black medv
## alto 325.7682 20.33720
## bajo 390.5477 25.31162
##
## Coefficients of linear discriminants:
## LD1
## zn 0.005331755
## nox -8.025497066
## age -0.016986785
## dis -0.053038787
## rad -0.074033448
## tax 0.001113865
## ptratio -0.062819427
## black 0.001227959
## medv -0.033685626
boston_lda_class <- tibble(pred = predict(boston_lda, newdata = boston_test)$class,
real = boston_test$crimen)
boston_lda_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo LDA con age e indus")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.891
## 2 kap binary 0.778
Los resultados son muy similares a los del modelo de regresión logística siendo ligeramente peores, sin embargo, la gráfica construida por el LDA ofrece una interpretabilidad mayor a la obtenida en el literal 10.
Modelo KNN (K-Nearest Neighbors)
boston_train_scaled <- boston_train %>%
dplyr::select(-crim, -crimen, -indus, -chas, -lstat) %>%
scale()
boston_test_scaled <- boston_test %>%
dplyr::select(-crim, -crimen, -indus, -chas, -lstat) %>%
scale()
boston_knn <- knn(train = boston_train_scaled, test = boston_test_scaled, cl = boston_train$crimen, k = 5)
boston_knn_class <- tibble(pred = boston_knn,
real = boston_test$crimen)
boston_knn_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo KNN con todos las variables")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.931
## 2 kap binary 0.861
El clasificador \(KNN\) tiene una tasa de clasificación correcta del 91.1% utilizando \(k = 3\) y todas las variables, sin embargo, al utilizar \(k = 5\) y las mismas variables que la regresión logística se obtiene un 93.1% de predicción correcta.
Sección 8.4 - Ejercicio 9
El conjunto de datos OJ contiene 1070 compras donde los clientes compraron un jugo de naranja de Citrus Hill o de Maid Orange Juice y se guardaron un número de características de los clientes.
a) Conjuntos de entrenamiento y validación
Crear un conjunto de entrenamiento con 800 observaciones y 270 de prueba
oj_train <- as_tibble(OJ) %>%
rownames_to_column() %>%
sample_n(size = 800)
oj_test <- as_tibble(OJ) %>%
rownames_to_column() %>%
anti_join(oj_train, by = "rowname")
oj_train <- oj_train %>% dplyr::select(-rowname)
oj_test <- oj_test %>% dplyr::select(-rowname)
dim(oj_train)## [1] 800 18
## [1] 270 18
b) Árbol de decisión
##
## Classification tree:
## tree(formula = Purchase ~ ., data = oj_train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 9
## Residual mean deviance: 0.7393 = 584.8 / 791
## Misclassification error rate: 0.16 = 128 / 800
El árbol resultante tiene 10 nodos terminales, la variable más importante es LoyalCH con una ponderación de 63, con una tasa de clasificación incorrecta de 14.88% en entrenamiento.
c) Interepretar el resultado de uno de los nodos
Interpretemos el nodo 1: En este nodo ingresan 800 observaciones del conjunto de entrenamiento, hay 311 observaciones de la clase CH y 489 de la clase MM, se realiza una partición a la izquierda con 512 observaciones y la derecha con 285 observaciones bajo el criterio \(LoyalCH < 0.48285\).
d) Gráfica del árbol y resultados
Los nodos utilizan las variables CH y PriceDiff para realizar la clasificación, los nodos terminales indican la clasificación final.
e) Predicciones, matriz de confusión y clasificación correcta
oj_tree_class <- tibble(pred = predict(oj_tree, type="class", newdata = oj_test),
real = oj_test$Purchase)
oj_tree_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo CART para clasificación")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.815
## 2 kap binary 0.620
Se obtiene una tasa de clasificación correcta del 80%. En general la mayor cantidad de errores se cometió al clasificar las variables CH como MM, esto ocurrió 29 veces.
f) Aplicar cv.tree() para determinar el tamaño óptimo del árbol
## $size
## [1] 9 6 5 2 1
##
## $dev
## [1] 144 144 138 182 289
##
## $k
## [1] -Inf 0 5 12 127
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
g) Gráfica con el tamaño del árbol y la validación cruzada
h) Árbol con la menor tasa de error de validación cruzada incorreta
De acuerdo a la gráfica anterior, el árbol con tamaño 5 parece obtener el menor error de validación cruzada. Realicemos la poda de acuerdo a esta inforamción y grafiquemos el nuevo árbol.
i) Creación del árbol podado
j) Comparación de los errores de clasificación en entrenamiento
oj_prune_class_train <- tibble(pred = predict(oj_prune, type="class"),
real = oj_train$Purchase)
oj_prune_class_train %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión en entrenamiento para el modelo CART podado para clasificación")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.834
## 2 kap binary 0.632
Se obtiene una clasificación correcta del 84.5% para el modelo CART podado, esto se compara con el modelo completo que da una clasificación correcta del 85.2%, es decir, el modelo con validación cruzada es ligeramente peor en el conjunto de entrenamiento pero tiene una mayor capacidad de generalización gracias a que evita sobreajustar los datos.
oj_tree_class_train <- tibble(pred = predict(oj_tree, type="class"),
real = oj_train$Purchase)
oj_tree_class_train %>%
metrics(truth = real, estimate = pred)## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.84
## 2 kap binary 0.645
k) Comparación los resultados del árbol podado vs el completo en validación
oj_prune_class <- tibble(pred = predict(oj_prune, type="class", newdata = oj_test),
real = oj_test$Purchase)
oj_prune_class %>%
conf_mat(truth = real, estimate = pred) %>%
autoplot(type = "heatmap") +
labs(title = "Matriz de confusión para el modelo CART podado para clasificación")## # A tibble: 2 x 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.819
## 2 kap binary 0.628
En el árbol sin podar obtuvimos una tasa de clasificación correcta del 80%, sin embargo, al realizar la poda se obtuvo una tasa de clasificación correcta del 80.4%, que podría ser una mejora apreciable dependiendo del volumen de ventas de jugo de naranja asociado que trae la predicción correcta.